adv_vis_1

Чтение данных

data <- readRDS("data/raw/very_low_birthweight.RDS") 

# Добавление id 

data <- data %>%
  mutate(id = row_number())

# Перевод категориальных переменных в факторы

data <- data %>% mutate(
  across(c(id, race, inout, delivery, pvh, ivh, ipe, sex), ~ as.factor(.x))) 

Удаление колонок с более 100 пропусков и строк с пропусками

cleaned_data <- data %>% select(where(~ sum(is.na(.)) <= 100)) %>% 
  filter(across(everything(), ~ !is.na(.)))

## Warning: Using `across()` in `filter()` was deprecated in dplyr 1.0.8.
## ℹ Please use `if_any()` or `if_all()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Удаление выбросов

outliers <- function(x) {
  (x < quantile(x, 0.25) - 3*IQR(x)) | 
    (x > quantile(x, 0.75) + 3*IQR(x))
}

cleaned_data %>% filter(across(where(is.numeric), ~ !outliers(.))) %>% glimpse()

## Warning: Using `across()` in `filter()` was deprecated in dplyr 1.0.8.
## ℹ Please use `if_any()` or `if_all()` instead.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Rows: 277
## Columns: 20
## $ birth    <dbl> 81.514, 81.558, 81.610, 81.624, 81.626, 81.689, 81.697, 81.70…
## $ exit     <dbl> 81.539, 81.667, 81.697, 81.700, 81.730, 81.878, 81.952, 81.82…
## $ hospstay <int> 9, 40, 32, 28, 38, 69, 93, 44, 44, 70, 85, 58, 75, 35, 34, 46…
## $ lowph    <dbl> 7.250000, 7.250000, 7.320000, 7.160000, 7.039997, 7.419998, 7…
## $ pltct    <int> 244, 182, 282, 153, 229, 361, 255, 186, 134, 229, 68, 174, 17…
## $ race     <fct> white, black, black, black, white, white, black, white, white…
## $ bwt      <int> 1370, 1480, 1255, 1350, 1310, 1180, 770, 1490, 1000, 1120, 74…
## $ gest     <dbl> 32.0, 32.0, 29.5, 34.0, 32.0, 28.0, 26.0, 33.0, 28.0, 29.0, 2…
## $ inout    <fct> born at Duke, born at Duke, born at Duke, born at Duke, born …
## $ twn      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ delivery <fct> abdominal, vaginal, vaginal, abdominal, vaginal, abdominal, v…
## $ apg1     <int> 7, 8, 9, 4, 6, 6, 4, 8, 5, 9, 9, 9, 1, 8, 4, 6, 9, 6, 8, 9, 8…
## $ vent     <int> 0, 0, 0, 0, 1, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1…
## $ pneumo   <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ pda      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld      <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0…
## $ year     <dbl> 81.51471, 81.55847, 81.61053, 81.62421, 81.62695, 81.68988, 8…
## $ sex      <fct> female, male, female, female, male, male, male, male, female,…
## $ dead     <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ id       <fct> 2, 4, 7, 10, 11, 14, 16, 17, 21, 22, 23, 25, 31, 35, 36, 41, …

Создание графиков плотности для числовых переменных

data_long <- cleaned_data %>%
  select(where(is.numeric)) %>%
  pivot_longer(everything(), names_to = "variable", values_to = "value")

ggplot(data_long, aes(x = value)) +
  geom_density(fill = "blue") +
  facet_wrap(~ variable, scales = "free")+
  theme_bw()

# Графики плотности для переменных bwt и gest с окраской по inout

data_long <- cleaned_data %>%
  pivot_longer(cols = c(bwt, gest), names_to = "variable", values_to = "value")

ggplot(data_long, aes(x = value, fill = inout)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ variable, scales = "free", nrow=2)+
  theme_bw()

# Cравнение значений колонки ‘lowph’ между группами в переменной inout

H0: среднее lowph в группе born at Duke равно среднему в группе transported

H1: средние не равны

Большие выборки, распределение приближается к нормальному по ЦПТ, используем двусторонний t тест c поправкой Уэлча для неизвестных дисперсий (уровень значимости установим 0.05)

t = 5.5731, df = 111.03, p-value = 1.77e-07

p-value<0.05

95 percent confidence interval: 0.05762847 0.1212193

0 не входит в ДИ

Следовательно, отклоняем нулевую гипотезу в пользу альтернативной

Как бы вы интерпретировали результат, если бы знали, что более низкое значение lowph ассоциировано с более низкой выживаемостью?

В таком случае группа transported имеет более низкую выживаемость по сравнению с группой born at Duke

test <- cleaned_data %>% t_test(lowph ~ inout, detailed=TRUE) %>% 
  add_significance() %>%
  mutate(y.position = max(cleaned_data$lowph, na.rm = TRUE) + 0.5) 

## Warning: The statistic is based on a difference or ratio; by default, for
## difference-based statistics, the explanatory variable is subtracted in the
## order "born at Duke" - "transported", or divided in the order "born at Duke" /
## "transported" for ratio-based statistics. To specify this order yourself,
## supply `order = c("born at Duke", "transported")`.

print(test)

## # A tibble: 1 × 8
##   statistic  t_df     p_value alternative estimate lower_ci upper_ci y.position
##       <dbl> <dbl>       <dbl> <chr>          <dbl>    <dbl>    <dbl>      <dbl>
## 1      5.57  111. 0.000000177 two.sided     0.0894   0.0576    0.121       8.05

  ggplot(cleaned_data, aes(x=inout, y=lowph, fill=inout))+
  geom_boxplot()+
  #stat_pvalue_manual(test, label="p-value")+ не работает с этой строкой, в чем причина?
  theme_bw()

Корреляционный анализ

num_data <- cleaned_data %>% select(where(is.numeric), -c(birth,year, exit)) %>% glimpse()

## Rows: 531
## Columns: 12
## $ hospstay <int> 9, 40, 2, 32, 28, 38, 62, 69, 1, 93, 44, 66, 65, 44, 70, 85, …
## $ lowph    <dbl> 7.250000, 7.250000, 6.969997, 7.320000, 7.160000, 7.039997, 7…
## $ pltct    <int> 244, 182, 54, 282, 153, 229, 182, 361, 378, 255, 186, 260, 18…
## $ bwt      <int> 1370, 1480, 925, 1255, 1350, 1310, 1110, 1180, 970, 770, 1490…
## $ gest     <dbl> 32.0, 32.0, 28.0, 29.5, 34.0, 32.0, 28.0, 28.0, 28.0, 26.0, 3…
## $ twn      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1…
## $ apg1     <int> 7, 8, 5, 9, 4, 6, 6, 6, 2, 4, 8, 1, 8, 5, 9, 9, 9, 6, 2, 1, 6…
## $ vent     <int> 0, 0, 1, 0, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 1, 0, 0, 1, 0, 1, 0…
## $ pneumo   <int> 0, 0, 1, 0, 0, 0, 1, 0, 1, 0, 0, 1, 1, 0, 0, 0, 0, 1, 0, 1, 0…
## $ pda      <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0…
## $ cld      <int> 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0…
## $ dead     <int> 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0…

data_cor <- cor(num_data)

data_cor %>% corr.test(method = "spearman") %>% 
  print(short=FALSE)

## Call:corr.test(x = ., method = "spearman")
## Correlation matrix 
##          hospstay lowph pltct   bwt  gest   twn  apg1  vent pneumo   pda   cld
## hospstay     1.00 -0.76 -0.71 -0.87 -0.83 -0.61 -0.76  0.84   0.64  0.83  0.95
## lowph       -0.76  1.00  0.83  0.85  0.90  0.60  0.84 -0.97  -0.94 -0.85 -0.85
## pltct       -0.71  0.83  1.00  0.78  0.75  0.41  0.91 -0.85  -0.83 -0.82 -0.80
## bwt         -0.87  0.85  0.78  1.00  0.97  0.80  0.92 -0.90  -0.74 -0.92 -0.95
## gest        -0.83  0.90  0.75  0.97  1.00  0.83  0.87 -0.94  -0.80 -0.94 -0.94
## twn         -0.61  0.60  0.41  0.80  0.83  1.00  0.65 -0.64  -0.39 -0.65 -0.71
## apg1        -0.76  0.84  0.91  0.92  0.87  0.65  1.00 -0.86  -0.76 -0.83 -0.86
## vent         0.84 -0.97 -0.85 -0.90 -0.94 -0.64 -0.86  1.00   0.88  0.91  0.93
## pneumo       0.64 -0.94 -0.83 -0.74 -0.80 -0.39 -0.76  0.88   1.00  0.83  0.73
## pda          0.83 -0.85 -0.82 -0.92 -0.94 -0.65 -0.83  0.91   0.83  1.00  0.94
## cld          0.95 -0.85 -0.80 -0.95 -0.94 -0.71 -0.86  0.93   0.73  0.94  1.00
## dead         0.69 -0.97 -0.79 -0.82 -0.87 -0.55 -0.80  0.90   0.97  0.85  0.78
##           dead
## hospstay  0.69
## lowph    -0.97
## pltct    -0.79
## bwt      -0.82
## gest     -0.87
## twn      -0.55
## apg1     -0.80
## vent      0.90
## pneumo    0.97
## pda       0.85
## cld       0.78
## dead      1.00
## Sample Size 
## [1] 12
## Probability values (Entries above the diagonal are adjusted for multiple tests.) 
##          hospstay lowph pltct  bwt gest  twn apg1 vent pneumo  pda  cld dead
## hospstay     0.00  0.07  0.11 0.01 0.03 0.20 0.08 0.02   0.20 0.03 0.00 0.13
## lowph        0.00  0.00  0.03 0.02 0.00 0.20 0.02 0.00   0.00 0.02 0.02 0.00
## pltct        0.01  0.00  0.00 0.06 0.08 0.37 0.00 0.02   0.03 0.03 0.05 0.05
## bwt          0.00  0.00  0.00 0.00 0.00 0.05 0.00 0.00   0.08 0.00 0.00 0.03
## gest         0.00  0.00  0.01 0.00 0.00 0.03 0.01 0.00   0.05 0.00 0.00 0.01
## twn          0.04  0.04  0.18 0.00 0.00 0.00 0.20 0.20   0.37 0.20 0.11 0.20
## apg1         0.00  0.00  0.00 0.00 0.00 0.02 0.00 0.01   0.08 0.03 0.01 0.05
## vent         0.00  0.00  0.00 0.00 0.00 0.02 0.00 0.00   0.01 0.00 0.00 0.00
## pneumo       0.03  0.00  0.00 0.01 0.00 0.21 0.00 0.00   0.00 0.03 0.10 0.00
## pda          0.00  0.00  0.00 0.00 0.00 0.02 0.00 0.00   0.00 0.00 0.00 0.02
## cld          0.00  0.00  0.00 0.00 0.00 0.01 0.00 0.00   0.01 0.00 0.00 0.06
## dead         0.01  0.00  0.00 0.00 0.00 0.06 0.00 0.00   0.00 0.00 0.00 0.00
## 
##  Confidence intervals based upon normal theory.  To get bootstrapped values, try cor.ci
##             raw.lower raw.r raw.upper raw.p lower.adj upper.adj
## hspst-lowph     -0.93 -0.76     -0.33  0.00     -0.96      0.00
## hspst-pltct     -0.91 -0.71     -0.22  0.01     -0.95      0.07
## hspst-bwt       -0.96 -0.87     -0.58  0.00     -0.98     -0.23
## hspst-gest      -0.95 -0.83     -0.48  0.00     -0.98     -0.13
## hspst-twn       -0.88 -0.61     -0.05  0.04     -0.92      0.15
## hspst-apg1      -0.93 -0.76     -0.32  0.00     -0.96      0.00
## hspst-vent       0.51  0.84      0.95  0.00      0.15      0.98
## hspst-pneum      0.10  0.64      0.89  0.03     -0.13      0.93
## hspst-pda        0.49  0.83      0.95  0.00      0.14      0.98
## hspst-cld        0.83  0.95      0.99  0.00      0.62      0.99
## hspst-dead       0.20  0.69      0.91  0.01     -0.08      0.95
## lowph-pltct      0.49  0.83      0.95  0.00      0.14      0.98
## lowph-bwt        0.55  0.85      0.96  0.00      0.19      0.98
## lowph-gest       0.68  0.90      0.97  0.00      0.37      0.99
## lowph-twn        0.04  0.60      0.87  0.04     -0.14      0.91
## lowph-apg1       0.51  0.84      0.95  0.00      0.15      0.98
## lowph-vent      -0.99 -0.97     -0.88  0.00     -1.00     -0.71
## lowph-pneum     -0.98 -0.94     -0.81  0.00     -0.99     -0.58
## lowph-pda       -0.96 -0.85     -0.55  0.00     -0.98     -0.19
## lowph-cld       -0.96 -0.85     -0.53  0.00     -0.98     -0.17
## lowph-dead      -0.99 -0.97     -0.88  0.00     -1.00     -0.71
## pltct-bwt        0.36  0.78      0.93  0.00      0.03      0.97
## pltct-gest       0.31  0.75      0.92  0.01     -0.01      0.96
## pltct-twn       -0.21  0.41      0.80  0.18     -0.30      0.83
## pltct-apg1       0.70  0.91      0.97  0.00      0.40      0.99
## pltct-vent      -0.96 -0.85     -0.55  0.00     -0.98     -0.19
## pltct-pneum     -0.95 -0.83     -0.48  0.00     -0.98     -0.13
## pltct-pda       -0.95 -0.82     -0.46  0.00     -0.98     -0.12
## pltct-cld       -0.94 -0.80     -0.41  0.00     -0.97     -0.07
## pltct-dead      -0.94 -0.79     -0.40  0.00     -0.97     -0.06
## bwt-gest         0.90  0.97      0.99  0.00      0.76      1.00
## bwt-twn          0.41  0.80      0.94  0.00      0.07      0.97
## bwt-apg1         0.74  0.92      0.98  0.00      0.47      0.99
## bwt-vent        -0.97 -0.90     -0.66  0.00     -0.99     -0.34
## bwt-pneum       -0.92 -0.74     -0.29  0.01     -0.96      0.02
## bwt-pda         -0.98 -0.92     -0.74  0.00     -0.99     -0.47
## bwt-cld         -0.99 -0.95     -0.83  0.00     -0.99     -0.62
## bwt-dead        -0.95 -0.82     -0.46  0.00     -0.98     -0.11
## gest-twn         0.48  0.83      0.95  0.00      0.12      0.98
## gest-apg1        0.58  0.87      0.96  0.00      0.23      0.98
## gest-vent       -0.98 -0.94     -0.79  0.00     -0.99     -0.54
## gest-pneum      -0.94 -0.80     -0.41  0.00     -0.97     -0.06
## gest-pda        -0.98 -0.94     -0.79  0.00     -0.99     -0.54
## gest-cld        -0.98 -0.94     -0.79  0.00     -0.99     -0.54
## gest-dead       -0.96 -0.87     -0.58  0.00     -0.98     -0.23
## twn-apg1         0.12  0.65      0.89  0.02     -0.13      0.93
## twn-vent        -0.89 -0.64     -0.11  0.02     -0.93      0.13
## twn-pneum       -0.79 -0.39      0.24  0.21     -0.79      0.24
## twn-pda         -0.89 -0.65     -0.12  0.02     -0.94      0.15
## twn-cld         -0.91 -0.71     -0.24  0.01     -0.95      0.06
## twn-dead        -0.86 -0.55      0.03  0.06     -0.89      0.17
## apg1-vent       -0.96 -0.86     -0.57  0.00     -0.98     -0.21
## apg1-pneum      -0.93 -0.76     -0.32  0.00     -0.96      0.01
## apg1-pda        -0.95 -0.83     -0.49  0.00     -0.98     -0.13
## apg1-cld        -0.96 -0.86     -0.57  0.00     -0.98     -0.21
## apg1-dead       -0.94 -0.80     -0.41  0.00     -0.97     -0.06
## vent-pneum       0.62  0.88      0.97  0.00      0.28      0.99
## vent-pda         0.70  0.91      0.97  0.00      0.40      0.99
## vent-cld         0.76  0.93      0.98  0.00      0.50      0.99
## vent-dead        0.66  0.90      0.97  0.00      0.34      0.99
## pneum-pda        0.48  0.83      0.95  0.00      0.12      0.98
## pneum-cld        0.26  0.73      0.92  0.01     -0.04      0.96
## pneum-dead       0.88  0.97      0.99  0.00      0.71      1.00
## pda-cld          0.79  0.94      0.98  0.00      0.54      0.99
## pda-dead         0.53  0.85      0.96  0.00      0.17      0.98
## cld-dead         0.36  0.78      0.93  0.00      0.03      0.97

corrplot(data_cor, method = 'number',insig = "blank")

data_cor %>%
  network_plot(min_cor = .0)

# Иерархическая кластеризация

# Стандартизация
data_scaled <- scale(data_cor)
# Матрица дистанций
data_dist <- dist(data_scaled, method = "euclidean")
# Дендрограмма кластеров
data_dist.hc <- hclust(d = data_dist, method = "ward.D2")
# Визуализация
fviz_dend(data_dist.hc, 
          k=2, 
          cex = 0.6,
          color_labels_by_k = TRUE,
          rect = TRUE) 

## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## ℹ The deprecated feature was likely used in the factoextra package.
##   Please report the issue at <https://github.com/kassambara/factoextra/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Определение оптимального количества кластеров методом силуэт и потом оптимизация кода выше
fviz_nbclust(data_scaled, kmeans, method = "silhouette")

# Heatmap+кластеризация

Вывод: в данных имеется два кластера: 1(pneumo, dead, hospstay, pda, vent, cld) и 2 (lowph, bwt, gest, pltct, apg1, twn). Между кластерами наблюдается преимущественно отрицательная корреляция.

pheatmap(data_scaled, 
         show_rownames = TRUE, 
         clustering_distance_rows = data_dist,
         clustering_method = "ward.D2", 
         cutree_rows = 5,
         cutree_cols = length(colnames(data_scaled)),
         angle_col = 45, 
         main = "Dendrograms for clustering rows and columns with heatmap")

# PCA

Вывод: первая компонента описывает 65% вариаций, а первые две компоненты описывают 75.45% вариаций данных (следовательно, в этих данных много скоррелированных между собой переменных). 1 группа (pneumo, dead, hospstay, pda, vent, cld) и 2 (lowph, bwt, gest, pltct, apg1, twn) имеют отрицательную корреляцию между группами, но при этом они скоррелированы между собой внутри групп.

# Используем стандартизованные значения для PCA, так как шкалы у переменных разные

data_pca <- prcomp(num_data, scale = T) 
summary(data_pca)

## Importance of components:
##                           PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.9667 1.1625 1.07954 0.99768 0.90725 0.88700 0.83410
## Proportion of Variance 0.3223 0.1126 0.09712 0.08295 0.06859 0.06556 0.05798
## Cumulative Proportion  0.3223 0.4350 0.53208 0.61503 0.68362 0.74918 0.80716
##                            PC8     PC9    PC10    PC11    PC12
## Standard deviation     0.81596 0.76417 0.66649 0.60125 0.50858
## Proportion of Variance 0.05548 0.04866 0.03702 0.03013 0.02155
## Cumulative Proportion  0.86264 0.91130 0.94832 0.97845 1.00000

fviz_pca_var(data_pca, col.var = "contrib")

# biplots

biplot <- ggbiplot(data_pca,
         groups = as.factor(cleaned_data$dead),
         scale=0, alpha = 0.5) + 
  theme_bw()
biplot

data_pca$id <-cleaned_data$id

ggplotly(biplot) %>%
  style(text = paste("ID:", data_pca$id)) # как сделать без перекрытия переменных?

Колонку dead использовать некорректно, так как она не учитывает время жизни (exit-birth), что необходимо для модели выживаемости

UMAP

Вывод: Два кластера (голубой и красный в углу) отличаются от остальных по каким-то признакам, при этом часть живых пациентов имеют характеристики, близкие к умершим. То же самое видно на PCA (умершие имеют позитивную корреляцию с одним кластером переменных и обратную корреляцию с другим кластером переменных, но в то же время часть живых пациентов тоже имеют такую зависимость).

umap_prep <- recipe(~., data = num_data) %>% 
  step_normalize(all_predictors()) %>% 
  step_umap(all_predictors()) %>%  
  prep() %>%   
  juice() 

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + 
  geom_point(aes(color = as.character(num_data$dead)),
             alpha = 0.7, size = 1) +
  labs(color = NULL) 

# Изменим n_neighbors и min_dist

Вывод: появилось 2 голубых кластера, данные выглядят более однородными

umap_prep <- recipe(~., data = num_data) %>% 
  step_normalize(all_predictors()) %>% 
  step_umap(all_predictors(),neighbors=5, min_dist=0.5) %>%  # по умолчанию было 15 и 0.01 соответственно
  prep() %>%   
  juice() 

umap_prep %>%
  ggplot(aes(UMAP1, UMAP2)) + 
  geom_point(aes(color = as.character(num_data$dead)),
             alpha = 0.7, size = 1) +
  labs(color = NULL)